{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "# 1.5 Spaces and forms on subdomains\n",
    "     \n",
    "\n",
    "In NGSolve, finite element spaces can be defined on subdomains. This is useful for multiphysics problems like fluid-structure interaction.\n",
    "\n",
    "In addition, bilinear or linear forms can be defined as integrals over regions. *Regions* are parts of the domain. They may be subdomains, or  parts of the domain boundary, or parts of subdomain interfaces.\n",
    "\n",
    "In this tutorial, you will learn about \n",
    "\n",
    "- defining finite element spaces on `Region`s,\n",
    "- `Compress`ing such spaces,\n",
    "- integrating on `Regions`s, and \n",
    "- set operations on `Region` objects via `Mask`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "from netgen.geom2d import *\n",
    "from ngsolve import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Naming subdomains and boundaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We define a geometry with multiple regions and assign names to these regions below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "geo = SplineGeometry()\n",
    "geo.AddRectangle((0,0), (2,2),\n",
    "                 bcs=[\"b\",\"r\",\"t\",\"l\"],\n",
    "                 leftdomain=1)\n",
    "geo.AddRectangle((1,0.9), (1.3,1.4),\n",
    "                 bcs=[\"b2\",\"r2\",\"t2\",\"l2\"],\n",
    "                 leftdomain=2, rightdomain=1)\n",
    "geo.SetMaterial (1, \"outer\")\n",
    "geo.SetMaterial (2, \"inner\")\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "\n",
    "* These statements define two rectangular subdomain regions named \"inner\", \"outer\". (Note how the ordering in `SetMaterial` is 1-based, not 0-based.)\n",
    "\n",
    "* The bottom, right, top and left parts of the outer rectangle's boundaries define boundary regions, respectively labeled  \"b\", \"r\", \"t\", \"l\".\n",
    "\n",
    "* Similarly, the bottom, right, top and left parts of the inner rectangle's boundaries are regions named \"b2\", \"r2\", \"t2\", \"l2\".\n",
    "\n",
    "* When you select the `Mesh` tab in the Netgen window's Graphical User Interface (GUI) and double click on a point, the two subdomains are rendered in different colors. (You should also see the subdomain outlines in the Netgen window when you select the `Geometry` tab.)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### A finite element space on a subdomain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes1 = H1(mesh, definedon=\"inner\")\n",
    "u1 = GridFunction(fes1, \"u1\")\n",
    "u1.Set(x*y)\n",
    "Draw(u1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Note how $u_1$ is displayed only in the `inner` region in the Netgen GUI."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Do not be confused with `ndof` for spaces on subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = H1(mesh)\n",
    "fes1.ndof, fes.ndof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Although `ndof`s are the same for spaces on one  subdomain and the whole domain, `FreeDofs` show that the real number of degrees of freedom for the subdomain space is much smaller:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(fes1.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "One can also glean this information by examining `CouplingType` of each degrees of freedom of the space, which reveals that there are many unknowns of type `COUPLING_TYPE.UNUSED_DOF`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(fes1.ndof):\n",
    "    print(fes1.CouplingType(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "It is possible to remove these dofs (thus making `GridFunction` vectors smaller) as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fescomp = Compress(fes1)\n",
    "print(fescomp.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Integrating on regions "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "You have already seen how boundary regions are used in setting Dirichlet boundary conditions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"b|l|r\")\n",
    "\n",
    "u, v = fes.TnT()\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Such boundary regions or subdomains can also serve as domains of integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(fes) \n",
    "f += u1 * v * dx(definedon=mesh.Materials(\"inner\")) \n",
    "f += 0.1 * v * ds(definedon=mesh.Boundaries(\"t\"))\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Here the functional $f$ is defined as a sum of two integrals, one over the inner subdomain, and another over the top boundary:\n",
    "$$\n",
    "f(v) = \\int_{\\Omega_{inner}} u_1 v \\; dx + \\int_{\\Gamma_{top}} \\frac{v}{10}\\; ds\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The python objects `dx` and `ds` represent volume and boundary integration, respectively. They admit a keyword argument `definedon` that is **either** a `Region` **or** a `str`, as can be seen from the documentation: look for \n",
    "`definedon: Optional[Union[ngsolve.comp.Region, str]]` in the help output below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "help(dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Thus \n",
    "- `dx(definedon=mesh.Materials(\"inner\"))` may be replaced by \n",
    "`dx('inner')`,\n",
    "\n",
    "and similarly, \n",
    "\n",
    "- `ds(definedon=mesh.Boundaries(\"t\")` may be replaced by  `ds('t')`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = LinearForm(fes) \n",
    "f += u1*v*dx('inner') + 0.1*v*ds('t')\n",
    "f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx\n",
    "a.Assemble()\n",
    "\n",
    "# Solve the problem:\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### More about `Region` objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.GetMaterials()    # list all subdomains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.GetBoundaries()   # list boundary/interface regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.Materials(\"inner\") # look at object's type "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "mesh.Boundaries(\"t\")    # same type!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Operations with `Region`s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Print region information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(mesh.Materials(\"inner\").Mask())\n",
    "print(mesh.Materials(\"[a-z]*\").Mask())  # can use regexp\n",
    "print(mesh.Boundaries('t|l').Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Add regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "io = mesh.Materials(\"inner\") + mesh.Materials(\"outer\")\n",
    "print(io.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Take complement of a region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "c = ~mesh.Materials(\"inner\")\n",
    "print(c.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Subtract regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "diff = mesh.Materials(\"inner|outer\") - mesh.Materials(\"outer\")\n",
    "print(diff.Mask())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Set a piecewise constant `CoefficientFunction` using the subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "domain_values = {'inner': 3.7,  'outer': 1}\n",
    "values_list = [domain_values[mat]\n",
    "               for mat in mesh.GetMaterials()]\n",
    "cf = CoefficientFunction(values_list)\n",
    "Draw(cf, mesh, 'piecewise')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Coefficients on boundary regions are given similarly using\n",
    "`mesh.GetBoundaries()`.\n",
    "\n",
    "### Extra exercises to try out\n",
    "\n",
    "Let's make a linear function that equals 2 at the bottom right vertex of current domain (0,2) x (0,2) and equals zero at the remaining vertices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "bdry_values = {'b': x, 'r': 2-y}\n",
    "values_list = [bdry_values[bc]\n",
    "               if bc in bdry_values.keys() else 0\n",
    "               for bc in mesh.GetBoundaries()]\n",
    "cf = CoefficientFunction(values_list)\n",
    "Draw(cf, mesh, 'piecewise')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "##### Pitfall!\n",
    "\n",
    "Look at the GUI output. This is not the function we wanted!\n",
    "\n",
    "What happened here? The object `cf` has no information on boundary regions, so it cannot associate the list of values to boundary regions. By default, the list of values are assumed to be subdomain values.\n",
    "\n",
    "To associate these values to boundary regions, we use a `GridFunction` and `Set`, which also lets us view an extension of these boundary values into the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g = GridFunction(H1(mesh), name='bdry')\n",
    "g.Set(cf, definedon=~mesh.Boundaries(''))\n",
    "Draw(g, max=2, min=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "If you think that specifying the whole boundary using\n",
    "\n",
    "    ~mesh.Boundaries('')\n",
    "\n",
    "is convoluted, then you may use\n",
    "\n",
    "    g.Set(cf, definedon=mesh.Boundaries('b|r')).\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g.Set(cf, definedon=mesh.Boundaries('b|r'))\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "What happens if you `Set` using `b` and then on `r`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g.Set(cf, definedon=mesh.Boundaries('b'))\n",
    "g.Set(cf, definedon=mesh.Boundaries('r'))\n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "##### Pitfall!\n",
    "Do you see an unexpected behavior in the GUI?\n",
    "\n",
    "If so, realize that each `Set` begins by zeroing the grid function - so any previous values are lost!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "What happens if you use\n",
    "\n",
    "    g.Set(cf, BND)\n",
    "    \n",
    "instead?\n",
    "\n",
    "##### Pitfall!   \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g = GridFunction(H1(mesh))\n",
    "g.Set(cf, BND)\n",
    "Draw(g, max=2, min=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "We do not obtain the previous result because the keyword `BND` directs `Set` to only set the data in boundary regions marked as `dirichlet`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.hycell": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "g = GridFunction(H1(mesh, dirichlet='b|r'))\n",
    "g.Set(cf, BND)\n",
    "Draw(g, max=2, min=0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "name": "subdomains.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
